<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>zcgesv.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="ZCGESV.1"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a>(N,NRHS,A,LDA,IPIV,B,LDB,X,LDX,WORK,SWORK,ITER,
     +                  INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK PROTOTYPE driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. WARNING: PROTOTYPE ..
</span><span class="comment">*</span><span class="comment">     This is an LAPACK PROTOTYPE routine which means that the
</span><span class="comment">*</span><span class="comment">     interface of this routine is likely to be changed in the future
</span><span class="comment">*</span><span class="comment">     based on community feedback.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER INFO,ITER,LDA,LDB,LDX,N,NRHS
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER IPIV(*)
      COMPLEX SWORK(*)
      COMPLEX*16 A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*)
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="ZCGESV.27"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a> computes the solution to a real system of linear equations
</span><span class="comment">*</span><span class="comment">     A * X = B,
</span><span class="comment">*</span><span class="comment">  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="ZCGESV.31"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a> first attempts to factorize the matrix in SINGLE COMPLEX PRECISION
</span><span class="comment">*</span><span class="comment">  and use this factorization within an iterative refinement procedure to
</span><span class="comment">*</span><span class="comment">  produce a solution with DOUBLE COMPLEX PRECISION normwise backward error
</span><span class="comment">*</span><span class="comment">  quality (see below). If the approach fails the method switches to a
</span><span class="comment">*</span><span class="comment">  DOUBLE COMPLEX PRECISION factorization and solve.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The iterative refinement is not going to be a winning strategy if
</span><span class="comment">*</span><span class="comment">  the ratio SINGLE PRECISION performance over DOUBLE PRECISION performance
</span><span class="comment">*</span><span class="comment">  is too small. A reasonable strategy should take the number of right-hand
</span><span class="comment">*</span><span class="comment">  sides and the size of the matrix into account. This might be done with a 
</span><span class="comment">*</span><span class="comment">  call to <a name="ILAENV.41"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> in the future. Up to now, we always try iterative refinement.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The iterative refinement process is stopped if
</span><span class="comment">*</span><span class="comment">      ITER &gt; ITERMAX
</span><span class="comment">*</span><span class="comment">  or for all the RHS we have:
</span><span class="comment">*</span><span class="comment">      RNRM &lt; SQRT(N)*XNRM*ANRM*EPS*BWDMAX 
</span><span class="comment">*</span><span class="comment">  where
</span><span class="comment">*</span><span class="comment">      o ITER is the number of the current iteration in the iterative
</span><span class="comment">*</span><span class="comment">        refinement process
</span><span class="comment">*</span><span class="comment">      o RNRM is the infinity-norm of the residual
</span><span class="comment">*</span><span class="comment">      o XNRM is the infinity-norm of the solution
</span><span class="comment">*</span><span class="comment">      o ANRM is the infinity-operator-norm of the matrix A
</span><span class="comment">*</span><span class="comment">      o EPS is the machine epsilon returned by <a name="DLAMCH.53"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>('Epsilon')
</span><span class="comment">*</span><span class="comment">  The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of linear equations, i.e., the order of the
</span><span class="comment">*</span><span class="comment">          matrix A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NRHS    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of right hand sides, i.e., the number of columns
</span><span class="comment">*</span><span class="comment">          of the matrix B.  NRHS &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input or input/ouptut) COMPLEX*16 array,
</span><span class="comment">*</span><span class="comment">          dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the N-by-N coefficient matrix A.
</span><span class="comment">*</span><span class="comment">          On exit, if iterative refinement has been successfully used
</span><span class="comment">*</span><span class="comment">          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
</span><span class="comment">*</span><span class="comment">          unchanged, if double precision factorization has been used
</span><span class="comment">*</span><span class="comment">          (INFO.EQ.0 and ITER.LT.0, see description below), then the
</span><span class="comment">*</span><span class="comment">          array A contains the factors L and U from the factorization
</span><span class="comment">*</span><span class="comment">          A = P*L*U; the unit diagonal elements of L are not stored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IPIV    (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The pivot indices that define the permutation matrix P;
</span><span class="comment">*</span><span class="comment">          row i of the matrix was interchanged with row IPIV(i).
</span><span class="comment">*</span><span class="comment">          Corresponds either to the single precision factorization 
</span><span class="comment">*</span><span class="comment">          (if INFO.EQ.0 and ITER.GE.0) or the double precision 
</span><span class="comment">*</span><span class="comment">          factorization (if INFO.EQ.0 and ITER.LT.0).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment">          The N-by-NRHS matrix of right hand side matrix B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B.  LDB &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
</span><span class="comment">*</span><span class="comment">          If INFO = 0, the N-by-NRHS solution matrix X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDX     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array X.  LDX &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) COMPLEX*16 array, dimension (N*NRHS)
</span><span class="comment">*</span><span class="comment">          This array is used to hold the residual vectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SWORK   (workspace) COMPLEX array, dimension (N*(N+NRHS))
</span><span class="comment">*</span><span class="comment">          This array is used to use the single precision matrix and the 
</span><span class="comment">*</span><span class="comment">          right-hand sides or solutions in single precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ITER    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          &lt; 0: iterative refinement has failed, double precision
</span><span class="comment">*</span><span class="comment">               factorization has been performed
</span><span class="comment">*</span><span class="comment">               -1 : taking into account machine parameters, N, NRHS, it
</span><span class="comment">*</span><span class="comment">                    is a priori not worth working in SINGLE PRECISION
</span><span class="comment">*</span><span class="comment">               -2 : overflow of an entry when moving from double to
</span><span class="comment">*</span><span class="comment">                    SINGLE PRECISION
</span><span class="comment">*</span><span class="comment">               -3 : failure of <a name="SGETRF.113"></a><a href="sgetrf.f.html#SGETRF.1">SGETRF</a>
</span><span class="comment">*</span><span class="comment">               -31: stop the iterative refinement after the 30th
</span><span class="comment">*</span><span class="comment">                    iterations
</span><span class="comment">*</span><span class="comment">          &gt; 0: iterative refinement has been sucessfully used.
</span><span class="comment">*</span><span class="comment">               Returns the number of iterations
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  if INFO = i, U(i,i) computed in DOUBLE PRECISION is
</span><span class="comment">*</span><span class="comment">                exactly zero.  The factorization has been completed,
</span><span class="comment">*</span><span class="comment">                but the factor U is exactly singular, so the solution
</span><span class="comment">*</span><span class="comment">                could not be computed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      COMPLEX*16 NEGONE,ONE
      PARAMETER (NEGONE=(-1.0D+00,0.0D+00),ONE=(1.0D+00,0.0D+00))
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL DOITREF
      INTEGER I,IITER,ITERMAX,OK,PTSA,PTSX
      DOUBLE PRECISION ANRM,BWDMAX,CTE,EPS,RNRM,XNRM
      COMPLEX*16 ZDUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL <a name="CGETRS.140"></a><a href="cgetrs.f.html#CGETRS.1">CGETRS</a>,<a name="CGETRF.140"></a><a href="cgetrf.f.html#CGETRF.1">CGETRF</a>,<a name="CLAG2Z.140"></a><a href="clag2z.f.html#CLAG2Z.1">CLAG2Z</a>,<a name="XERBLA.140"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>,ZAXPY,
     $         ZGEMM,<a name="ZLACPY.141"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>,<a name="ZLAG2C.141"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER IZAMAX
      DOUBLE PRECISION <a name="DLAMCH.145"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>,<a name="ZLANGE.145"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
      EXTERNAL IZAMAX,<a name="DLAMCH.146"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>,<a name="ZLANGE.146"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC ABS,DBLE,MAX,SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      DOUBLE PRECISION   CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      ITERMAX = 30
      BWDMAX = 1.0E+00
      DOITREF = .TRUE.
<span class="comment">*</span><span class="comment">
</span>      OK = 0
      INFO = 0
      ITER = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      IF (N.LT.0) THEN
          INFO = -1
      ELSE IF (NRHS.LT.0) THEN
          INFO = -2
      ELSE IF (LDA.LT.MAX(1,N)) THEN
          INFO = -4
      ELSE IF (LDB.LT.MAX(1,N)) THEN
          INFO = -7
      ELSE IF (LDX.LT.MAX(1,N)) THEN
          INFO = -9
      END IF
      IF (INFO.NE.0) THEN
          CALL <a name="XERBLA.181"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>(<span class="string">'<a name="ZCGESV.181"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a>'</span>,-INFO)
          RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if (N.EQ.0).
</span><span class="comment">*</span><span class="comment">
</span>      IF (N.EQ.0) RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Skip single precision iterative refinement if a priori slower
</span><span class="comment">*</span><span class="comment">     than double precision factorization.
</span><span class="comment">*</span><span class="comment">
</span>      IF (.NOT.DOITREF) THEN
          ITER = -1
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute some constants.
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="ZLANGE.199"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>(<span class="string">'I'</span>,N,N,A,LDA,WORK)
      EPS = <a name="DLAMCH.200"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>(<span class="string">'Epsilon'</span>)
      CTE = ANRM*EPS*SQRT(DBLE(N))*BWDMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set the pointers PTSA, PTSX for referencing SA and SX in SWORK.
</span><span class="comment">*</span><span class="comment">
</span>      PTSA = 1
      PTSX = PTSA + N*N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert B from double precision to single precision and store the
</span><span class="comment">*</span><span class="comment">     result in SX.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLAG2C.211"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,NRHS,B,LDB,SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -2
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert A from double precision to single precision and store the
</span><span class="comment">*</span><span class="comment">     result in SA.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLAG2C.221"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,N,A,LDA,SWORK(PTSA),N,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -2
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the LU factorization of SA.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGETRF.230"></a><a href="cgetrf.f.html#CGETRF.1">CGETRF</a>(N,N,SWORK(PTSA),N,IPIV,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -3
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve the system SA*SX = SB.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGETRS.239"></a><a href="cgetrs.f.html#CGETRS.1">CGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,SWORK(PTSA),N,IPIV,
     +            SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert SX back to double precision
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CLAG2Z.244"></a><a href="clag2z.f.html#CLAG2Z.1">CLAG2Z</a>(N,NRHS,SWORK(PTSX),N,X,LDX,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute R = B - AX (R is WORK).
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLACPY.248"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>      CALL ZGEMM(<span class="string">'No Transpose'</span>,<span class="string">'No Transpose'</span>,N,NRHS,N,NEGONE,A,LDA,X,
     +           LDX,ONE,WORK,N)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check whether the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">     stopping criterion. If yes, set ITER=0 and return.
</span><span class="comment">*</span><span class="comment">
</span>      DO I = 1,NRHS
          XNRM = CABS1(X(IZAMAX(N,X(1,I),1),I))
          RNRM = CABS1(WORK(IZAMAX(N,WORK(1,I),1),I))
          IF (RNRM.GT.XNRM*CTE) GOTO 10
      END DO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If we are here, the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">     stopping criterion. We are good to exit.
</span><span class="comment">*</span><span class="comment">
</span>      ITER = 0
      RETURN
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 30 IITER = 1,ITERMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Convert R (in WORK) from double precision to single precision
</span><span class="comment">*</span><span class="comment">         and store the result in SX.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="ZLAG2C.275"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,NRHS,WORK,N,SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span>          IF (INFO.NE.0) THEN
              ITER = -2
              GO TO 40
          END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Solve the system SA*SX = SR.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="CGETRS.284"></a><a href="cgetrs.f.html#CGETRS.1">CGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,SWORK(PTSA),N,IPIV,
     +                SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Convert SX back to double precision and update the current
</span><span class="comment">*</span><span class="comment">         iterate.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="CLAG2Z.290"></a><a href="clag2z.f.html#CLAG2Z.1">CLAG2Z</a>(N,NRHS,SWORK(PTSX),N,WORK,N,INFO)
<span class="comment">*</span><span class="comment">
</span>          CALL ZAXPY(N*NRHS,ONE,WORK,1,X,1)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Compute R = B - AX (R is WORK).
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="ZLACPY.296"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>          CALL ZGEMM(<span class="string">'No Transpose'</span>,<span class="string">'No Transpose'</span>,N,NRHS,N,NEGONE,A,
     +               LDA,X,LDX,ONE,WORK,N)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Check whether the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">         stopping criterion. If yes, set ITER=IITER&gt;0 and return.
</span><span class="comment">*</span><span class="comment">
</span>          DO I = 1,NRHS
              XNRM = CABS1(X(IZAMAX(N,X(1,I),1),I))
              RNRM = CABS1(WORK(IZAMAX(N,WORK(1,I),1),I))
              IF (RNRM.GT.XNRM*CTE) GOTO 20
          END DO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         If we are here, the NRHS normwised backward errors satisfy the 
</span><span class="comment">*</span><span class="comment">         stopping criterion, we are good to exit.
</span><span class="comment">*</span><span class="comment">
</span>          ITER = IITER
<span class="comment">*</span><span class="comment">
</span>          RETURN
<span class="comment">*</span><span class="comment">
</span>   20     CONTINUE
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If we are at this place of the code, this is because we have
</span><span class="comment">*</span><span class="comment">     performed ITER=ITERMAX iterations and never satisified the stopping
</span><span class="comment">*</span><span class="comment">     criterion, set up the ITER flag accordingly and follow up on double
</span><span class="comment">*</span><span class="comment">     precision routine.
</span><span class="comment">*</span><span class="comment">
</span>      ITER = -ITERMAX - 1
<span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Single-precision iterative refinement failed to converge to a
</span><span class="comment">*</span><span class="comment">     satisfactory solution, so we resort to double precision.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZGETRF.333"></a><a href="zgetrf.f.html#ZGETRF.1">ZGETRF</a>(N,N,A,LDA,IPIV,INFO)
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLACPY.335"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,X,LDX)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.EQ.0) THEN
          CALL <a name="ZGETRS.338"></a><a href="zgetrs.f.html#ZGETRS.1">ZGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,A,LDA,IPIV,X,LDX,INFO)
      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZCGESV.343"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a>.
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
